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1—5 , Abstract 

00 . 

Special high-accuracy direct force summation A^-body algorithms and their rele- 

, ^ vance for the simulation of the dynamical evolution of star clusters and other gravi- 

.^ ' tating A^-body systems in astrophysics are presented, explained and compared with 

lO . other methods. Other methods means here approximate physical models based on 

'^ \ the Fokker-Planck equation as well as other, approximate algorithms to compute 

f— ^ ■ the gravitational potential in A'^-body systems. Questions regarding the parallel im- 

0> . plementation of direct "brute force" A^-body codes are discussed. The astrophysical 

application of the models to the theory of relaxing rotating and non-rotating colli- 

sional star clusters is presented, briefly mentioning the questions of the validity of 

the Fokker-Planck approximation, the existence of gravothermal oscillations and of 



O \ rotation and primordial binaries. 



/\ ' 1 Introduction 



"The dynamical evolution of an isolated spherical system composed of very 
many mass points has an appealing simplicity. The Newtonian laws of motion 
are exact, and all average quantities are functions only of radial distance r 
and time t. Nevertheless, it is only recently, with the availability of fast com- 
puters, that a systematic understanding of how such systems develop through 
time has emerged. Since these idealized systems should provide a very good 
approximation for globular clusters in this and other galaxies, the theory of 
their development is an important part of astronomy as well as an interesting 
branch of theoretical particle dynamics." [77] 



* to appear in: Riffert H., Werner K. (eds), Computational Astrophysics, The 
Journal of Computational and Applied Mathematics (JCAM), Elsevier Press, 
Amsterdam. 



Preprint submitted to Elsevier Preprint 1 February 2008 



Once celestial mechanics was one of the most important fields of astronomy. 
Nowadays astrophysics has become much wider in scope, including fields like 
stellar astrophysics and gas and plasma dynamics of interstellar matter. For 
some objects, however, the pure dynamics of gravitating mass points still 
provides an excellent description of the global dynamical evolution or gives 
at least the dominating background in which the gaseous or baryonic matter 
evolves. Such objects are, starting from the large scale, the entire universe 
itself, some evolutionary phases of galaxies and galactic nuclei, globular and 
open star clusters, and last but not least our planetary system. Globular star 
clusters are gas free systems with some 10^ stars, orbiting around our own [14] 
and other galaxies [70]. 

This article aims at the complex interplay of thermodynamic processes like 
heat conduction and relaxation with the physics of self-gravitating systems 
and the stochastic nature of star clusters having finite particle number A^, 
and the specific computational and physical models used for the numerical 
simulation of the dynamical evolution of star clusters under these processes 
on the computer. Globular clusters are a very good laboratory for relaxation 
processes in discrete particle systems, because their dynamical and relaxation 
timescales are well separated from each other and from the lifetime of the clus- 
ter and of the universe as a whole. In this article the methods appropriate to 
model their evolution are in the focus. Other kinds of iV-body simulations are 
useful for example for hydrodynamics ("smoothed particle hydrodynamics"), 
galaxy dynamics ( "collisionless systems" ) or cosmological A^-body simulations 
of structure formation in the universe and are covered by other articles in this 
volume. The main distinction of those from the models presented here, is 
that the dynamics of systems dominated by two-body relaxation ( "collisional 
systems") requires typically very high accuracies (typical energy error per 
crossing time AE/ E < 10^^ or smaller) over very long physical integration 
times (thousands of crossing times). The term "collisional" here always refers 
to systems, whose evolution is infiuenced by relaxation through elastic two- or 
more-body encounters, not to physical collisions, where two stars collide and 
merge or disrupt each other. As a consequence of the high accuracy require- 
ments for collisional A^-body simulations, commonly known algorithms like 
the leap-frog time integration and the TREE-method to compute the gravi- 
tational potential of a particle distribution, are not efficient to use here; the 
use of high-order time integration schemes and "brute-force" algorithms to 
compute the potential are more efficient, as will be argued below. 

This article is organized as follows: this introduction is followed by a section on 
the approximate models of self-gravitating collisional A^-body systems, after 
which practical and theoretical aspects of the corresponding highly accurate 
direct A^-body simulations are presented. Finally astrophysical applications of 
the methods and relevant questions under study are presented. 



Let us begin with the definition of some useful time scales. A typical particle 
crossing time tcr in a star cluster is 



f = — 



where Vh is the radius containing 50 % of the total mass and ah is a typical 
velocity associated with the root mean square random motion (velocity disper- 
sion) taken at r/j. If virial equilibrium prevails, we have a^ ~ GMh/vh (where 
the sign ^ here and henceforth means "approximately equal" or "equal within 
an order of magnitude"), thus 



GMh 



(2) 



This is equal to the dynamical timescale, which is also used for example in 
the theory of stellar structure and evolution. Global dynamical adjustments 
of the system, like oscillations, are connected with this timescale. Taking the 
square of equation 2 yields t^j, ~ r\/{GMh) which is related to Kepler's third 
law, because the orbital velocity in a Keplerian point mass potential has the 
same order of magnitude as the velocity dispersion in virial equilibrium. 

Unlike most laboratory gases stellar systems are not usually in thermodynamic 
equilibrium, neither locally nor globally. Radii of stars are usually extremely 
small relative to the average interparticle distances of stellar systems (e.g. the 
radius of the sun is r© ^ 10^° cm, a typical distance between stars in our 
galactic neighbourhood is of the order of lO^^cm). Only under rather special 
conditions in the centres of galactic nuclei and during the short high-density 
core collapse phase of a globular cluster, stellar densities might become large 
enough that stars come close enough to each other to collide, merge or disrupt 
each other. 

Therefore it is extremely unlikely under normal conditions that two stars 
touch each other during an encounter; encounters or collisions usually are 
elastic gravitative scatterings. Fairly generally the mean interparticle distance 
is large compared to po = Gm/o"^, which is the impact parameter for a 90° 
defiection in a typical encounter of two stars of equal mass m, where the rel- 
ative velocity at infinity is \/2(T, with local ID velocity dispersion a. Thus 
most encounters are small-angle defiections. The relaxation time trx is defined 
as the time after which the root mean square velocity increment due to such 
small angle gravitative defiections is of the same order as the initial velocity 
dispersion of the system. We use the local relaxation time as defined by [49]: 

_ 9 a^ 



G is the gravitational constant, p the mean stellar mass density, A^ the to- 
tal particle number, and 7 a parameter of order unity, which results from 
an integration over all possible impact parameters for a two-body encounter. 
Taking the linear system dimension as a maximum impact parameter yields 
7 = 0.4 [77]. Measurements in direct star by star evolutionary simulations of 
stellar systems are more in favour of 7 = 0.11 [22,23], which is the value used 
throughout this article. 

Assuming virial equilibrium a fundamental proportionality turns out: 

0^ 1 / AA • (4) 
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(cf. e.g. [77]). As a result, for very large A^, dynamical equilibrium is attained 
much faster than thermodynamic equilibrium. If one assumes a purely kinetic 
temperature definition, it ensues that in star clusters the temperatures (or 
velocity dispersions) can remain different for different coordinate directions 
over many dynamical times. For example, in a spherical system the radial 
and tangential velocity dispersion would be different, which is denoted as 
anisotropy. 

There are several reasons to believe that anisotropy is present and important 
for the dynamical evolution of astrophysical star clusters. Many observations 
are matched better by models including anisotropy [55], and all direct sim- 
ulations exhibit the formation of anisotropy under very general conditions, 
independent of the underlying physical cause driving the system's evolution. 
[8] showed in the context of a gas dynamical model of star clusters, that 
isotropy remains only under very special conditions (linear profiles of veloc- 
ities of mass and energy transport), and a similar study [36] gives the same 
result for axisymmetric collapsing gaseous systems. 



2 Approximate Models 

2. 1 Fokker-Planck Approximation 



Unfortunately, the direct simulation of such rich stellar systems as globular 
clusters with star-by-star modelling is not yet possible. The gap between the 
largest useful A^-body models with particle numbers of the order of a few 10*^ 
particles and the median globular star cluster (A^ ~ 5 x 10^) can only be 
bridged at present by use of theory. There are two main classes of theory: 
(i) Fokker-Planck models, which are based on the Boltzmann equation of the 



kinetic theory of gases [13,72,29,15], and (ii) gas models [54,79,32], which can 
be thought of as a set of moment equations of the Fokker-Planck model. 

These simplified models are the only detailed models which are directly appli- 
cable to large systems such as globular clusters. But their simplicity stems from 
many approximations and assumptions which are required in their formula- 
tion. Examples are the assumptions of spherical symmetry, which contradicts 
the asymmetry of the galactic tidal field, or statistical estimates of cross sec- 
tions for the formation of close binaries by three-body or dissipative (tidal) 
two-body encounters, and for their subsequent gravitational interactions with 
field stars. Such processes play a dominant role to reverse core collapse of 
globular clusters, which otherwise would inevitably lead to a singular density 
profile with infinite density at the centre [9,19,41]. 

The Fokker-Planck approximation truncates the so-called B^GKY hierarchy 
of kinetic equations (see [10]) at lowest order assuming that for most of the 
time all particles are uncorrelated with each other and only coupled via the 
smooth global gravitational potential. Correlations only play a role as a se- 
quence of uncorrelated two-body encounters. Instead of determining a general 
correlation function one resorts to a phenomenological description of the ef- 
fects of collisions by computing diffusion coefficients directly from the known 
solution of the two-body problems. Diffusion coefficients D{Avi) and D{AviVj) 
denote the average rate of change of Vi and ViVj due to the cumulative effect of 
many small angle deffections during two-body encounters. Let m, v and rrif, 
Vf be the mass and velocity of a star from a test and field star distribution, re- 
spectively (both distributions can but need not to be the same). In Cartesian 
geometry the diffusion coefficients are defined by 

D{Avi)=47iG^mf\nA^^ ; D{Av,v,) = AnG^rriflnA^j^ ; (5) 



where / is the phase space density of stars (briefiy: distribution function) and 
g, h are the Rosenbluth potentials defined in [76] 

h{v) = {m + raf) / — — -^d^Vf ; g{v) = rrif / f {vf)\v — Vf\(Pvf . (6) 

i/ (J ^ f \ 



Note that provided the distribution function / is given in terms of a convenient 
polynomial series as in Legendre polynomials the Rosenbluth potentials can 
be evaluated analytically to arbitrary order, as was seen already by [76], see 
for a modern rederivation and its use for star cluster dynamics [26,82]. With 
these results we can finally write down the local Fokker-Planck equation in its 
standard form for the Cartesian coordinate system of the vf. 
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The subscript "enc" should refer to encounters, which are the driving force 
of two-body relaxation. Still Eq. 7 is a six- dimensional integro-differential 
equation; its direct numerical simulation in stellar dynamics can presently 
only be done by further simplification. First Jeans' theorem is applied and / 
transformed into a function of the classical integrals of motion of a particle 
in a potential under the given symmetry, as e.g. energy E and modulus of 
the angular momentum J^ in a spherical potential or E and 2;-component of 
angular momentum J^ in axisymmetric coordinates. Thereafter the Fokker- 
Planck equation can be integrated over the accessible coordinate space for 
any given combination of constants of motion and the orbit-averaged Fokker- 
Planck equation ensues. By transformation from Vi to E and J and via the 
limits of the orbital integral the potential enters both implicitly and explicitly. 
In a two-step scheme alternatively solving the Poisson- and Fokker-Planck 
equation a direct numerical solution is obtained [12,15,85-87,18]. One of the 
main uncertainties in this method is that for non-spherical mass distributions 
the orbit structure in the system may depend on unknown non-classical third 
integrals of motion which are neglected. 



2.2 Anisotropic Gaseous Model 



The local Fokker-Planck equation Eq. 7 is utilized in another way for gaseous 
or conducting sphere models of star clusters. Integrating it over velocity space 
with varying powers of the velocity coordinates yields a system of equations in 
the spatial coordinates; the local approximation is used in the sense that the 
orbit structure of the system is not taken into account, diffusion coefficients 
and all other quantities are assumed to be well defined just as a function of 
the local quantities (density, velocity dispersions and so on). The system of 
moment equations is truncated in third order by a phenomenological equation 
of heat transfer. Such approach has been suggested by [56,33] and generalized 
to anisotropic systems by [8,54], and for a presentation of the recent model see 
e.g. [26]. In the following the derivation of the model equations is decribed. 



2.2.1 The "Left Hand Stdes" 

In spherical symmetry, polar coordinates r 6, (p are used and t denotes the 
time. The vector v = {vi),i = r, 9, 0, denotes the velocity in a local Cartesian 



coordinate system at the spatial point r,6,(f). For brevity u = Vr, v = vg, 
w = v^ is used. The distribution function /, which due to spherical symmetry 
is a function of r, t, u, v"^ + w"^ only, is normalized according to 



p{r,t) = / f{r,u,v + w ,t)dudvdw, (9) 



where p{r, t) is the mass density; if m denotes the stellar mass, we get the 
particle density n = p/m. Then 

u= uf{r,u,v'^ + w'^,t)dudvdw, (10) 



is the bulk radial velocity of the stars. Note that for the analogously defined 
quantities v and w we have v = w = 0. 

In order to go ahead to the anisotropic gaseous model equations we now turn 
back to the left hand side of the Fokker-Planck equation Eq. 7, which is the 
collisionless Boltzmann or Vlasov operator. For practical reasons we prefer for 
the left hand side local Cartesian velocity coordinates, whose axes are oriented 
towards the r, 6, coordinate space directions. With the Lagrange function 

the Euler-Lagrange equations of motion for a star moving in the cluster po- 
tential $ become: 

(9$ v^ + w^ . uv w"^ . uw vw ,_„, 



The complete local Fokker-Planck equation, derived from Eq. 7, attains the 
form 

^ + u^ + u^ + v^ + w^=(^J-] (13) 

dt dr du dv dw \ St I 

\ / one 



where the term subscribed by "enc" denotes the terms involving diffusion 
coefficients as in Eq. 8. Moments {i,j, k) of / are defined in the following way 
(all integrations range from — oo to oo): 



(0,0,0):=p= fdudvdw ; (1,0,0) :=«= ufdudvdw (14) 

{2,0,0): = pr + pu^ = u^fdudvdw (15) 
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(0,2,0) ■. = pg = v'^fdudvdw ; {0,0,2) := p^ = w^fdudvdw (16) 

(3, 0, 0) := Fr + 3upr + u^ = I u^fdudvdw (17) 

(1,2,0) := Fg + upe = fuv^fdudvdw (18) 

(1,0,2) := F^ + up^ = uw'^fdudvdw . (19) 



Note that the definitions of pi and Fi are such that they are proportional to 
the random motion of the stars. Due to spherical symmetry we have pg = 
p^=: Pt and Fg = F^ =: Ft/2. By Pr = pa^ and pt = pa^ the random velocity 
dispersions are given, which are closely related to observables in globular star 
clusters and galaxies. It is convenient to define velocities of energy transport 
by 

Fr Ft 

Vr = ^ + u ; vt = \-u . (20) 

3pr 2pt 



By multiplication of the Fokker-Planck equation 13 with various powers of u, 
V, w we get up to second order the following set of moment equations (for a 
detailed derivation in the here used variables see [78] , bar for u dropped in the 
following) : 



I + ^^IP"")-" (21) 

^ + „^ + £^ + i|£ + 2^^ = (22) 

at or r^ p or pr 

^ + 1 A (r^upr) + 2pr— + -— (r^Fr) - ^ = i^-^\ (23) 

dt r^ dr ^ '^' ^ dr r"^ Or ^ ^' r \ St I , . „ 

\ / cnc,bin3 

— + ——(r^u ]+2— + - — —(r^F]+ — =(—\ (24) 

dt r"^ dr ^ ' r 2r'^ dr ^ ' r \ 6t I , . „ 

\ / cnc,bin3 

The terms labeled with "enc" and "bin3" symbolically denote the collisional 
terms resulting from the moments of the right hand side of the Fokker-Planck 
equation (Eq. 8) and an energy generation by formation and hardening of 
three body encounters. Both will be discussed below. With the definition of 
the mass Mr contained in a sphere of radius r 

-—S. = Anry (25) 

or 



the set of Eqs. 22-24 is equivalent to gasdynamical equations coupled with 
Poisson's equation. Since moment equations of order n contain moments of 



order n+1, it is necessary to close the system of the above equations by 
an independent closure relation. Here we choose the heat conduction closure, 
which consists of a phenomenological ansatz in analogy to gas dynamics. It 
was first used (restricted to isotropy) by [56]. It is assumed that heat transport 
is proportional to the temperature gradient, where we use for the temperature 
gradient an average velocity dispersion o'^ = (a"^ + 2af)/3 and assume Vr = 
Vt (this latter closure was first introduced by [8]). Therefore, the last two 
equations to close our model are 

Vr-U+ 5— =0 5 Vr = Vt. (26) 

AirGpt^^ or 



With Eqs. 22-24, 25, and 26 we have now seven equations for our seven de- 
pendent variables Mr, p, u, pr, Pt, Vr, Vt- 



2.2.2 Binary Heating 

It was already early realized that in a star cluster with single stars under high 
density conditions, one or more strongly bound binaries form, which could 
dominate the further evolution [35,4]. This is a contradiction to the basic 
assumption underlying the Fokker-Planck equation, that the only correlations 
in the system are those produced by a sequence of uncorrelated small-angle 
gravitational encounters. Nevertheless [9] introduced a phenomenological heat 
source into their gaseous model equations, in order to describe the input of 
random kinetic energy ("heat") to the cluster by formation and hardening 
of so-called three-body binaries. The ansatz for the functional form of the 
heating term has been clarified and more thoroughly discussed by [27,34]. They 
describe a simple estimate for the rate of formation of binaries by close three- 
body encounters of single stars; in subsequent superelastic scatterings between 
the formed binary and field stars the binary will on average become harder, 
provided its binding energy is large compared to the mean temperature of the 
surrounding single stars. Surplus kinetic energy taken from the gravitational 
binding energy of the binary members goes to the field star and thus provides 
a heating source for the core of the cluster. There is an upper limit of the 
binary binding energy given by the condition that the recoil on the binary 
in a typical superelastic scattering due to conservation of linear momentum 
in the process leads to escape of the binary. As a result each binary after its 
formation supplies a certain amount of energy by three-body encounters to the 
system until it escpes. The resulting heating term is (isotropic binary heating 
assumed) : 

2a™.v(^)% (|t) ^m (27) 




Here a simple estimate using gravitational focusing and the probability that 
three particles come together have been employed. Cb is a constant of propor- 
tionality which ies expected to have a value between 75 and 90 for an equal 
mass system; for more details see the above cited papers. 



2.2.3 The "Right Hand Sides" 

All right hand sides of the moment equations 22-24 are calculated by multi- 
plying the right hand side (the encounter term) of the Fokker-Planck equation 
as it occurs in Eq. 8 with the appropriate powers of m, v and w and integrating 
over velocity space. There is only one non-trivial encounter term to be deter- 
mined for the collisional decay of anisotropy. It is self-consistently computed 
by assuming a certain Legendre series evaluation for / up to second order (i.e. 
including anisotropy) in the Appendix of [26], the result being {pa = Vr — Pt)'- 

_ P^ _ 10 _ 9 (T^ 



ta defined in the above equation denotes the characteristic decay time of 
anisotropy; trx is equivalent to the standard two-body relaxation time. The 
particular factors applied to it originate unambigously from the Fokker-Planck 
collisional term evaluation with the assumption of a certain normalization and 
functional form of / by a Legendre series. The procedure can be thoroughly 
followed in [49]. For the above result terms quadratic in pa have been omitted. 
Comparisons with direct iV-body simulations suggest a more general ansatz 

^'^ (29) 





^ata 



and it is shown that Aa = 0.1 provides the best results [26]. Sect. 4 describes 
some examples how well the gaseous and Fokker-Planck models describe a 
star cluster's evolution as compared to a direct A^-body simulation. There is 
no other way to check the theoretical models on the Fokker-Planck equation, 
because the timescale for exponential instability and deterministic chaos to 
occur in a self-gravitating star cluster consisting of many stars of equal or at 
least similar mass is of the same order as a crossing time [28]. There is no 
analytical or semianalytical general solution of the A^-body problem available 
in that case for the unperturbed problem. In contrast to this in the case of 
solar system studies there is a semianalytic secular theory [50]) to be compared 
with the direct orbit integrations (see e.g. [51]). Here, for the star cluster case, 
we only can rely on the comparison of the numerical solutions obtained from 
different physical models, as there are direct A^-body integrations and models 
based on the Fokker-Planck approximation. 
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3 Direct A^-body Simulations — Methods and Algorithms 

3.1 Introduction - Density and Potential Computation 



To integrate the orbits of particles in time under their mutual gravitational 
interaction the total gravitational potential at each particle's position is re- 
quired. Poisson's equation in integral form gives the potential $ generated at 
a point in coordinate space r due to a smooth mass distribution p{r) 

$(f) = -G I T^^d^r' . (30) 

J \r — r\ 



There are two fundamentally different methods to define the density distribu- 
tion as a function of a given particle distribution. The first is based on a mesh 
in coordinate space; particles are sampled on the mesh and their mass divided 
by the cell volume, which provides a local density. This method called particle- 
mesh requires for good statistics a sufficient number of particles in each cell. 
There is no or very little intrinsic particle-particle relaxation with this method, 
but there is relaxation of particle energies due to the finite resolution of the 
mesh (see [38], and for a more recent cosmological application compare [44] 
and references therein). Refinements, by which particles are smeared out by 
low-order interpolation formulae (e.g. cloud-in-cell, CIC) or the acceleration 
is interpolated within the cells (e.g. Superbox [20]) are possible and reduce 
mesh relaxation. 

The second method is based on the particles itself. A kernel function W{x, h) 
is defined, normalized by / Wd^x = 1, where h describes a typical length scale 
over which the influence of a particle decays. Therewith a sampled density ps 
is defined in a mesh-free way as 

p^{r)= fw{f-f',h)p{f)d^f (31) 



where /i is a characteristic smoothing length. A discrete particle distribution 
is given by p(r) = J2^{^ ~ ^j) with N particles distributed at positions Vj. 
Hence we get a sampled density from the discrete particle distribution as 

N 

p^(f) = Y, mjW{f- r", h) (32) 

.7 = 1 



As an estimate for the density and other thermodynamic quantities this method 
is used by "smoothed particle hydrodynamics" simulations [71], using kernel 
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functions W with compact support, which means they are non-zero only in a 
bounded volume. Here, it is not intended to explain this further, the reader 
is referred to the literature and other papers of this volume. If one takes a S- 
function for the kernel as well and puts this into the integral Poisson equation 
(30) Newton's law turns out again: 



N 



$(r) 



-gj: 



rrii 



r — Ti 



(33) 



Certainly this could have been written down immediately, but the above de- 
scription makes it easier to understand the relation of the direct potential 
summation to the other methods. In the following, the prototype A^-body in- 
tegration method using the above "brute-force" method, in a specific form 
called the Hermite scheme [60], shall be described in some more detail. It is 
the most commonly used method in the field of globular cluster dynamics and 
other studies requiring very high accuracies. 



3.2 The Hermite Scheme 



Assume a set of A^ particles with positions rj(to) and velocities Vi{tQ) {i = 
1, . . . , A^) is given at time t = to, and let us look at a selected test particle 
at f = -To = ^(to) and v = Vq = vito). Note that here and in the following 
the index i for the test particle i and also occasionally the index indicating 
the time to will be dropped for brevity; sums over j are to be understood to 
include all j with j ^ i, since there should be no self-interaction. Accelerations 
ao and their time derivatives ao are calculated explicitly: 






Oo 



J2Gm^ 



Vj 3(V, ■ R,)R, 



R^ 



R] 



(34) 



where Rj :- 
predictions. 



fj, Vj := V — Vj, Rj := \Rj\, Vj := \Vj\. By low order 



Xp(t) = -{t - to)^do + -(t - tofaQ + (t - tQ)v + x 
^pit) = 2^t - hfo'O + (^ - ^o)ao + V , 



(35) 
(36) 



new positions and velocities for all particles at t > to are calculated and used 
to determine a new acceleration and its derivative directly according to Eq. 
34 at t = ti, denoted by ai and di. On the other hand di and di can also be 
obtained from a Taylor series using higher derivatives of a at t = to: 
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100 200 



500 1000 2000 



Fig. 1. The relative energy error as the function of the number of steps. A time-step 
criterion using differences between predicted and corrected values is used, differ- 
ent from Eq. 43. Dotted curves are for Hermite schemes, solid curves for Aarseth 
schemes. The stepnumber p denotes the order of the integrator. From [57]. 



1 1 

0,1 = -{t- t^fa^Q + -(t - tofaf' + (t - to)ao + ao , 
o 2 



1 



\2-(3) 



?(2) 



di = -(t - to) 4 + (^ - ^0)4 ^ + do • 



(37) 
(38) 



If di and di is known from direct summation (from Eq. 34 using the predicted 
positions and velocities) one can invert the equations above to determine the 
unknown higher order derivatives of the acceleration at t = to for the test 
particle: 



-d(2) 
2 

1. 
— ( 

6 



^g(3)^^«o-Qi 



do — di 2do + di 
it - to) 
do + di 



it-hf 



{t-t,f {t-t,y 



(39) 
(40) 



This is the Hermite interpolation, which finally allows to correct positions and 
velocities at ti to high order from 



x{t) = Xp{t) + ^(t - to)^dE,') + ^(t - to)^d(3) 
v{t) = v,it) + ^(t - to)^d[f) + ^(t - to)^(i^ 



(41) 
(42) 
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Taking the time derivative of Eq. 42 it turns out that the error in the force 
calculation for this method is 0{At^), as opposed to the widely used leap-frog 
schemes, which have a force error of 0{At^). Additional errors induced by 
approximate potential calculations (particle mesh or Tree) create potentially 
even larger errors than that. In Fig. 1, however, it is shown that the above 
Hermite method used for a real A^-body integration sustains an error of O(At^) 
for the entire calculation. Many persons in the world know as Aarseth scheme 
(in particular the code version NB0DY5 [1]) an integrator of the same order as 
the Hermite scheme, but using only accelerations on four time points instead 
of a and a on two time points. As is shown in [57], the Aarseth scheme is 
0{At'^) as well, but for the same number of time steps the absolute value 
of the energy error (not its slope) is clearly smaller in the Hermite scheme. 
This means that for a given energy error the Hermite scheme allows timesteps 
which are larger by some factor of order unity depending on the parameters 
of the system under study. The Hermite scheme has been commonly adopted 
during the past years, because it needs less memory, and allows slightly larger 
timesteps. More importantly, after the addition of a hierarchical (as opposed 
to individual) time step scheme it is well suited for parallelization on modern 
special and general purpose high performance computers [81]. The timestep 
scheme will be discussed now. 



3.3 Choice of Timesteps - Parallelization 
[1] provides an empirical timestep criterion 



At 



laP 



, ri^^ ^ ^ . 43 



The error is governed by the choice of f], which in most practical applications 
is taken to be r^ = 0.01 — 0.04. It is instructive to compare this with the inverse 
square of the curvature k, of the curve a{t) in coordinate space 

1 = i±ffi (44) 

K \a^ >\ 



Clearly under certain conditions the time step choice Eq. 43 becomes similar 
to choosing the timestep according to the curvature of the acceleration curve; 
since it was determined just empirically, however, it cannot generally be related 
to the curvature expression above. In [57] a different time step criterion has 
been suggested, which appears simpler and more straightforwardly defined, 
and couples the timestep to the difference between predicted and corrected 
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Fig. 2. The logarithm. of-X', -the-totaLcompu-ting-time required to advance the iV-body 
system for one crossing time plotted as a function of the particle number A^ for the 
equal time step scheme Tgq, the individual time step scheme Tind and the Ah- 
mad-Cohen neighour scheme with two levels of individual time steps Tu- The unit 
of computing time is the time required to calculate the force between a pair of 
particles. The system is assumed to be homogeneous. From [61]. 

coordinates. The standard Aarseth time step criterion Eq. 43 has been used 
in most A^-body simulations so far (but compare the discussion in [84]). 

Since the position of all field particles can be determined at any time by the 
low-order prediction Eq. 36, the time step of each particle (which determines 
the time at which the corrector Eq. 42 is applied) can be freely chosen accord- 
ing to the local requirements of the test particle; the additional error induced 
due to the use of only predicted data for the full N sums of Eq. 34 is negligibly 
small, for the benefit of not being forced to keep all particles in lockstep. Such 
an individual time step scheme is in particular for non-homogeneous systems 
very advantageous, as was quantitatively pointed out by [61]. Particles in the 
high density core of a star clusters need to be updated much more often than 
particles on orbits very far from the centre. They show that the gain in com- 
putational speed due to the individual time step scheme (as compared to a 
lockstep scheme where all particles share the minimum required time step) is 
of the order N'^^^ for homogeneous and A^^ for strongly spatially structured 
systems; we show their results as Figs. 2, 3. 

For the purpose of vectorization and parallelization it is better not to have 
the particles continuously distributed on a time axis. Consequently, [58] uses 
a hierarchical scheme, still on the basis of Eq. 43; but a change of the timestep 
is considered only if that equation yields a variation of At compared to the 
last step by more than a factor of 2 (increase or decrease). If this is the 
case a variation by 2 is applied only. Thus in model units all timesteps are 
selected from the set {2~*|i = 0,...imax} with k = imax determined by the 
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Fig. 3. As Fig. 2, but for the system with a power-law density distribution p oc r^^'^^. 
From [61]. 

condition that Atmin > 2^*™'"'= for the minimum timestep Atmin determined 
from Eq. 43. For core collapse simulations of star clusters of a few ten thousand 
particles imax goes up to about 20; empirically and theoretically [61] Atmin oc 
A^~^'^, so for large A^ imax becomes larger, however, on the other hand, how 
large imax grows for fixed A^ depends on the selected criteria for so-called KS 
regularisation of perturbed two-body motion (see below). The implementation 
of the block step scheme indeed uses an even stronger condition than the 
above described one, it is demanded that not only the time steps, but also 
the individual accumulated times of each particles are commensurate with the 
timestep itself. This ensures that for any particle i and any time Tj = tj + 6ti 
all particles with Stj < Sti have for their own time Tj = tj + 6tj = Ti, where 
the last equality is the non-trivial one. Such procedure is important for the 
parallelization of the algorithm. For example it has as a consequence that 
at the big time steps always huge groups of particles are due for correction, 
sometimes even all particles (at the largest steps). Such scheme allows an 
efficient parallelization of all operations necessary for calculation of a and d and 
for the update of particle positions and velocities (corrections). Special purpose 
computers have been built tailored to the Hermite codes, which are denoted 
as HARP ( "Hermite Accelerator Pipeline" ) boards and stem from the bigger 
GRAPE-family [83,63]. Such HARP-boards have been made available also 
at some places outside Japan, including "Astronomisches Rechen-Institut" 
Heidelberg (for an application see e.g. [88]). 

Another refinement of the Hermite or Aarseth "brute force" method is the 
two-time step scheme, denoted as neighbour or Ahmad-Cohen scheme [5]. For 
each particle a neighbour radius is defined, and a and a are computed due to 
neighbours and non-neighbours separately. Similar to the Hermite scheme the 
higher derivatives are computed separately for the neighbour force (irregular 
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Fig. 4. Theoretical speedup (neglecting communication) of regular force calculation 
as a function of processor number for varying particle number N . The dashed line is 
the ideal maximum speed up which could be reached on a given processor number. 

force) and non-neighbour force (regular force). Computing two timesteps, an 
irregular small Atirr and a regular large Atj-cg , from these two force components 
by Eq. 43 yields a timestep ratio of 7 := Atj-cg/Atirr being in a typical range 
of 5-20 for N of the order 10^ to 10'^. The reason is that the regular force has 
much less fluctuations than the irregular force. The Ahmad-Cohen neighbour 
scheme is implemented in a self-regulated way, where at each regular time-step 
a new neighbour list is determined using a given neighbour radius r^j for each 
particle. If the neighbour number found is larger than the prescribed optimal 
neighbour number, the neighbour radius is increased or vice versa. In [1,61] 
more complicated algorithms to adjust the neighbour radius are described. On 
the contrary to [61], who find an optimal neighbour number of iV„^opt oc A^^'^ 
we find that adopting a constant neighbour number of the order of 20 — 50 
is sufficient at least up to A^ = 50000. The reason is that by using special 
purpose machines or parallelization for parts of the code, an optimal neighbour 
number is not well defined, so the neighbour number can be selected according 
to accuracy and efficiency requirements [81]. After each regular time step the 
new neighbour list is communicated along with the new particle positions to all 
processors of the parallel machine, thus making it possible to do the irregular 
time step in parallel as well. 

Using a two-time step or neighbour scheme again increases the computational 
speed of the entire integration by a factor of at least proportional to N^l'^ 
[57]. Both the regular and irregular timesteps are arranged in the hierarchi- 
cal, commensurable way, and the total inherent parallelism in the resulting 
algorithm is depicted in Figs. 4, 5 from [81] for the irregular and the regular 
step. One can see that even for moderate particle numbers of 10*^ particles 
some 512 processors could be used efficiently. Sometimes there are only very 
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Fig. 5. As Fig. 4, for the irregular (neigbour) force calculation. 



few particles in the smallest steps to be integrated, which one might consider 
as being very prohibitive for parallelization. However, due to the large num- 
ber of medium and large size blocks this effect is negligible for the overall 
performance. It causes however, the saturation in the curves in Figs. 4 and 5 
which defines the limit for the number of processors useful for a given particle 
number A^. By using more and more processors in the parallel execution one 
finds that the asymptotic scaling of the "brute force" A^-body problem can be 
reduced effectively to an A^ scaling (Fig. 6). But in our present implementa- 
tion the parallelization is done only according to parallel sections (do loops) 
in the code; there is no domain decomposition (distributing particles on the 
processor). Thus at the end of any timesteps new results have to be broadcast 
to all other processing units. A systolic algorithm is used for that which scales 
linearly in communication time with the number of processors. It is interesting 
to note an approach suggested by molecular dynamicists to use a new kind 
of hyper-systolic communication algorithm, which scales only by the square 
root of the processor number [52,53]. Presently we think that hyper-systolic 
algorithms can efficiently be used only if the sum over all particles for the 
acceleration and its time derivative (Eq. 34) should be directly parallelized. 
The number of interprocessor communications A^comm for the hyper-systolic 
algorithm is of the order Ny/npE] on the other hand our algorithm, which we 
would like to call here "parallel group execution algorithm" [81], has a scaling 
-^comm oc N'^^^upE, bccause only subgroups of particles, whose size scales with 
N"^/^ have to be communicated across the processor network. In other words, 
asymptotically (above some critical particle number as a function of npE, the 
hyper-systolic algorithm should lose against the parallel group execution al- 
gorithm. However, these questions have not yet been examined in detail, for 
example what the critical N really are and which algorithm is more efficient 
for practically useful particle numbers of today. This is subject of present and 
future work. 




Fig. 6. CPU time needed for one A^-body time unit as a function of particle number 
N using NB0DY6++ on the CRAY T3E. The collection of data points includes 
runs with varying average neighbour number and processor/pipeline number, start- 
ing from 8 for low N up to 512 for the largest A^, which are not individually dis- 
criminated in the figure. 

If the two-body force between any pair of particles becomes dominant their 
(perturbed) relative motion is integrated in special regularized coordinates 
(taking into account perturbations from field particles), in which the singular- 
ity of the two-body motion is transformed into a slowly varying parameter (the 
binding energy) and does not occur in the integration variables. The rest of 
the A^-body simulation generally regards the regularized pair as a compound 
particle located at the position and moving with the velocity of its centre 
of mass, except in the case when a perturber moves very close to a regular- 
ized pair (in such cases the pair is resolved). It was already discovered in the 
earliest published A^-body simulations that the formation of close and eccen- 
tric binaries occurs as the rule rather than as an exception and that it was 
particularly difficult to accurately integrate them [39,40]. As a consequence 
two-body, three-body and chain regularizations were developed and imple- 
mented in order to accurately and efficiently integrate star clusters including 
all their close binaries, triples and hierarchical subsystems. An excellent ac- 
count of regularization, historically and scientifically, can be found in [66]. 
Most recent developments are the slow-down treatment of tight binaries [67] 
and a new method to gain accuracy and exact solutions in the unperturbed 
case using Stumpff functions [68] . 

Recently the necessity of regularization was challenged and its replacement 
by a binary tree structure for hierarchical systems with relative coordinates 
has been suggested [62]. However, the regularisation procedure is undisputedly 
much more efficient and accurate for highly eccentric binaries, and the new 
method has not yet been widely applied and proven to work through the 
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most delicate phase of core bounce and post-collapse evolution in point-mass 
systems or systems with many primordial hard binaries. 



3.4 Exponential Instability, Validity of N-body simulations, Planetary Sys- 
tem Integrations 



Concerning the algorithms explained in the previous paragraph the direct A^- 
body simulation may turn out to be the most reliable (although computation- 
ally most expensive) way to simulate the dynamical evolution of a gravitating 
system consisting of A^ point masses. It does not involve any serious approx- 
imations and assumptions, as e.g. the Fokker-Planck approximation in the 
gaseous models. By reducing the ?7-values any accuracy can be achieved in 
principle, as far as the globally conserved quantities (energy, angular momen- 
tum) are concerned. However, for a system with A^ particles phase space has 
6A^ dimensions, and a check of say energy and angular momentum alone only 
checks whether the numerically calculated system remains within the allowed 
6A^ — 4 dimensional hypervolume. There is no a priori information how "ex- 
act" the individual trajectories are reproduced in the simulation. [69] pointed 
out that, due to repeated close encounters occurring between particles initial 
configurations that are very close to each other, quickly diverge in their evo- 
lution from each other. He could show that the separation in phase space of 
two trajectories increases exponentially with time, or with other words, the 
evolution of the configuration is extremely sensitive to initial conditions (par- 
ticle positions and velocities). The timescale of exponential instability is as 
short as a fraction of a crossing time, and the accurate integration of a system 
to core collapse would require of order 0{N) decimal places [28,45]. Those 
papers argue that the problem is caused by two-body encounters, but chaotic 
orbits in non-integrable potentials can be a source of exponential instability 
and thus cause unreliable numerical integrations as well. 

However, the situation is not as bad as it seems. A^-body simulations for 
star clusters or galactic nuclei do not really exploit the detailed configuration 
space of all particles. Quantities of interest are global or somehow averaged 
quantities, like Lagrangian radii or velocity dispersions averaged in certain 
volumes. As it was nicely demonstrated in the pioneering series of papers by 
[22-25] such results are not sensitive to small variations of initial parameters. 
They took statistically independent initial models (positions and velocities at 
the beginning selected by different random number sets) and showed that the 
ensemble average of the dynamical evolution of the system always evolved 
predictably and in remarkable accord with results obtained from the Fokker- 
Planck approximation. The method was also partly and successfully used in 
[26] , which focused on the evolution of anisotropy and comparisons with the 
anisotropic gaseous models of the author of this paper. 



20 



As a consequence, it should be remembered, however, that great care has to 
be taken when interpreting results of A^-body simulations on a particle by 
particle basis, for example determining rates of specific types of encounters, 
which could produce mergers in a large direct iV-body model. 

The long-term behaviour of dynamical systems as the solar system are be- 
ing studied by iV-body simulations as well, but clearly there are much higher 
requirements on the accuracy of the individual orbits in contrast to the star 
cluster problem. Therefore for the solar system dynamics symplectic methods, 
using a generalized leap-frog, like the widely used Wisdom-Holman symplec- 
tic mapping method [89] are the standard integration method. As a non- 
exhaustive reference the reader might look into a recent study of the relation 
between the earth-moon system and the stability of the inner solar system 
using this method [43] and a contemporary review [16]. Symplectic mapping 
methods do not show secular errors in energy and angular momentum. How- 
ever, in their standard implementation they require a constant timestep. A 
generalization using a time transformation simultaneously with the general- 
ized leap-frog has been suggested which can cope with variable timesteps [65]. 
Another more practical approach to strongly reduce secular errors is to en- 
force a time-symmetric scheme by making the timesteps reversible through an 
iteration [42,21]. How well this generally works and its relation to symplectic 
schemes is presently not clear. In [68] it is stressed that even with a newly ap- 
plied classical method secular errors in the integration of close binaries can be 
strongly reduced. One should keep in mind though, that the A^-body integra- 
tion schemes discussed in this paper yield excellent results in the star cluster 
research (see Sect. 4) but are unsuitable for long-term solar system studies, 
because they generally have secular errors, although small. As outlined above 
in star cluster simulations the secular errors are being kept small relative to 
typical values of energy and angular momentum and an accurate reproduction 
of all individual stellar orbits is not generally required. 



3.5 What about TREE- and fast multipole codes? 



Finally, remarks shall be made on two very widely used algorithms to compute 
gravitational potentials from particle distributions namely the Tree- and fast 
multipole (FMP) algorithms. The Tree- method of [6] divides the system into 
hierarchical cells. The mutual interaction between particles or cells is resolved 
only if the opening parameter 9 = r/d, where r is the distance to and d a size 
scale of the cell under consideration, is smaller than a prescribed critical ^crit- 
If the cell is not resolved because 9 < ^crit, there is still the option to evaluate 
multipole moments of its internal mass distribution for the interaction with 
external particles. As one can see from Fig. 7, a global accuracy requirement of 
AE/E ^ 10"^ demands ^crit ~ 0.2, a value much smaller than the usually effi- 
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Fig. 7. Tradeoff between CPU time per step and average force error for the 
TREE-code with monopole terms only. From Fig. 4.11. of [74]. 
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Fig. 8. CPU time per step versus particle number N for TREE-codes with varying 
opening parameter 9 and a direct full "brute-force" labelled with PP in the figure 
(for "particle-particle"). From Fig. 4.9. of [74]. 

dent choice of 0.5-0.7, at which the computational time scales approximately 
as 0{N In N). Looking then at Fig. 8, the computational time for the TREE- 
code with ^crit ~ 0.2 scales nearly as 0{N'^), i.e. like a "brute force" algorithm. 
So for each particle number and required accuracy one should carefully check 
whether a TREE-code or a direct A^-body code are the best choice. 

Another TREE-based algorithm is the fast-multipole method (FMT) proposed 
by [30,31]. The pair- wise potential in Eq. 30 is approximated by a multipole 
series, which can be done for arbitrary precision if enough terms are included. 
The multipole terms used for different test particles can be transformed into 
each other by using clever addition theorems for spherical harmonics, so the 
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Fig. 9. Comparison of computational time as a function of particle number N be- 
tween particle-particle, hierarchical Tree, and the fast multipole code. From Fig. 
7.14 of [74]. 

entire algorithm scales in its computational demand with 0{N) only. Higher 
precision only changes the proportionality factor, not the scaling (as in the 
case of the TREE-code, which effectively becomes a "brute force" code if high 
enough accuracy is demanded. However, such a code is fine only for homoge- 
neous or nearly homogeneous systems, as they occur in plasma physics. In all 
cases where there is strong spatial structure, like in astrophysical star clus- 
ters, [61] have demonstrated that the use of an individual time step scheme 
in an 0{N'^) code gains a factor at least oc iV in efficiency. So, asymptotically 
a "brute-force" integrator with individual timesteps is more efficient than an 
FMT integrator. The latter is based on an equal timestep for all particles 
(otherwise it would lose its 0{N) property; so both codes have asymptoti- 
cally the same A^ scaling, but then the overhead (proportionality) factor is 
much smaller in the direct force summation than in the multipole evaluation. 
This can be seen also from Fig. 9 in [74] for low A^. For the direct calculation 
method in this plot, the individual time step scheme is not taken into account. 

The information contained in the previous paragraphs, complemented by some 
additional details and references, which will not be elaborated in more detail 
here, are presented in an overview in Table 1. It is divided into three boxes, 
the first for the mesh or series evaluation codes, which do not contain particle- 
particle forces and thus are not appropriate for direct modelling of relaxing 
systems. The second box contains the classical direct "brute force" A^-body 
codes, whereas the third one contains algorithms which cannot clearly be 
counted to one of the other two groups. 
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Table 1 



Algorithms for A^-body Simulations 

N: particle, Nn'- characteristic neighbour number 

ric'- number of grid cells in one dimension, nlm: order in 3D series evaluation 



Acronym 


Algorithm 


Scaling 


Comments 


PM 


Particle Mesh 


N nl log2 nl (1) 


fixed geometry 


FMP 


Fast Multipole 


N nlm 


req. equal At 


SCF 


Self-Consistent Field 


N nlm 


series evaluation (^' 


NbodyI 


Aarseth 


Ar2 


ITS, softening 


NBODY1 + + 


Hermite 


iV2 


HTS, softening 


Nbody2 


Aarseth, AC 


NNn + N^h 


ITS, softening, (3) 


Nbody3 


Aarseth 


N^ 


ITS, KS-reg. 


Nbody4 


Hermite 


iV2 


HTS, KS-reg. 


NbodyS 


Aarseth, AC 


NNn + A^Vt 


ITS, KS-reg., (^^ 


Nbody6 


Hermite, AC 


NNn + A^Vt 


HTS, KS-reg., (3) 


NBODY6 + + 


parallel Nbody6 


NN^ + Ny-f 


HTS, KS-reg., (3,4) 


KiRA 


Hermite 


N^ 


HTS, (5) 


Tree 


TREE-code 


NlnN 


A^^ for high accuracy 


F^M 


Part.-Part. PM 


A^2^3i„g^„3(l) 


fixed geometry ^^> 



softening: singularity in pairwise potential removed by softening parameter e 

ITS: Individual Time Step Scheme 

HTS: Hierarchical Block Time Step Scheme 

KS-reg.: KS regularization of perturbed two- and hierarchical A^-body motion [48,68] 

AC: Ahmad-Cohen neighbour scheme [5] 

^^' Discrete FFT on regular 3D mesh with n linear mesh points assumed 

(^) Sufficient Accuracy requires appropriate basis function set [37] 

(^^ 7: ratio of regular to irregular time step 

^^> speedup by parallel execution not contained in scaling, see [81] 

(^' New high accuracy Hermite code based on Starlab [64,75] 

'^^ with hierarchically nested adaptive grids used for cosmological simulations [73] 



4 Application to Star Clusters 



Since this article is focused on the physical and numerical methods of calculat- 
ing the evolution of relaxing star clusters, only a brief account of some of the 
physical problems and challenges will be given here, which have been and will 
be tackled by the previously described models. Despite of a wealth of beau- 
tiful observational data provided by e.g. Hubble space telescope observations 
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Fig. 10. Evolution of the 1% Lagrangian radius in the averaged N = 1000 iV-body 
model in comparison to the anisotropic gaseous model for different strength of the 
binary energy generation parameter Cb (see Eq. 27). The subscript cc indicates a 
pure core collapse gaseous model without binary heating. 
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Fig. 11. Evolution of the 1% Lagrangian radius in the averaged 1000, 50, 250 body 
models in comparison to the anisotropic gaseous models using C5 = 90, 70, 55, 
respectively. 

of globular clusters some of the fundamental questions related to the validity 
of the iV-body approach and the other approximate models still deserve at- 
tention as they can lead to very fascinating general questions regarding the 
thermodynamical behaviour of large A^-body systems. 

A series of papers has been devoted to the comparison of ensemble averaged 
A^-body simulations (A^ < 2000) with the expectations derived from Fokker- 
Planck or gaseous models [26,22-25]. Here we show as an example, in Figs. 
10 and 11, the excellent agreement reached between the anisotropic gaseous 
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Fig. 12. Lagrangian radii containing the indicated fraction of total mass as a func- 
tion of time for a single 10^-body simulation (fluctuating curves) compared to an 
averaged N = 1000 simulation of [22,23]. Times scaled as explained in main text. 

model and the ensemble averaged A^-body system. The models started with 
an initial Plummer model and follow the core collapse induced by heat con- 
duction and the post-collapse evolution due to formation and hardening of 
three-body hard binaries. The agreement of both types of models mutually 
supports both sides: it shows that by ensemble averaging, the exponential in- 
stability of the A^-body system does not spoil the physically correct behaviour 
of the system. It also demonstrates that the Fokker-Planck approximation, 
especially with its underlying assumption of strict spherical symmetry and 
dominance of small-angle two-body encounters for relaxation (i.e. neglectance 
of collective processes), is correct. It also shows that the very simple algorithm 
to describe the heating provided by the formation of close three-body binaries 
and their subsequent hardening by superelastic binary-single star encounters, 
which was first introduced into the gaseous models by [9], provides a surpris- 
ingly good description of the real processes in the average iV-body system. 
The cited paper ignited a discussion over many years whether gravothermal 
oscillations, being a thermodynamic consequence of heat conduction by two- 
body relaxation, will prevail in a real A^-body system with all its stochastic 
fluctuations. The question was settled after an A^-body simulation on the mas- 
sively parallel Teraflop GRAPE machine [63] using a high-accuracy Hermite 
scheme, as described above, became available. Gravothermal oscillations were 
found in a very large A^ = 32000 particle simulation [59] . 

In Fig. 12 we show a striking example of the validity of the Fokker-Planck 
approximation even for a single large direct A^-body-simulation, here using 
Nbody5, an Aarseth scheme (see Table 1), for 10000 particles, a model simu- 
lation again starting with Plummer's model and undergoing core collapse and 
core bounce due to hard binaries [80] . The average A^ = 1000 particle model by 
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Fig. 13. Evolution of central angular velocities versus central density of all 10 mass 
bins used in the calculation for an initial model with dimensionless central potential 
Wq = 3 and dimensionless angular velocity ujq = 0.30. The highest mass bin with 
m = 1.4 Mq is the uppermost curve, while the lowest mass bin with m = 0.15 Mq is 
the lowest curve. The parameters Wq and ujq refer to an initial Michie-King model. 

[22] has been taken and its time was scaled with the factor A^/ln(7A^), which 
is the scaling of the standard two-body relaxation time Eq. 4. An excellent 
match between the evolution of the Lagrangian radii for the two systems oc- 
curs after such scaling, proving that it is indeed the standard relaxation which 
dominates the pre-collapse evolution. The differences between the two systems 
show up at the moment of the formation of the first three-body binaries, after 
which one expects the evolution not to scale as the relaxation time. In simple 
terms, the larger A^, the less important are three-body effects as compared to 
the global potential; hence for large particle numbers the system collapses to 
higher densities and three-body effects finally dominate because they depend 
on the third power of the particle density as compared to the n^ dependence 
of two-body relaxation. 

Finally, we show a result from [17] in Fig. 13, a new multi-mass model using 
the orbit averaged Fokker-Planck approximation for axisymmetric rotating 
relaxing star clusters. The standard effect of mass segregation of the heavy 
masses is accompanied here by an acceleration of their rotational speed as 
compared to the small masses. Such interesting dynamical behaviour occurs 
just due to point-mass relaxation processes starting with a very simple tidally 
truncated rotating King model without any mass or rotational segregation. It 
is a yet unpublished generalization of equal mass rotating star cluster models 
[18]. They neglect the possible dynamical effect of non-classical third integrals, 
since it is assumed that the distribution function depends on energy and z- 
component of the angular momentum only. Such approximation needs to be 
checked by direct iV-body models, which is the subject of on-going work. The 
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results will also be important for the dynamical study of rotating galactic 
nuclei containing massive star-accreting black holes. 

The reader should be made aware of the problem of scaling in the description 
of escaping stars from globular clusters [75,3] being tackled by large direct 
A^-body simulations and their comparisons with approximate models. There 
are more challenges, like the inclusion of many close binaries already originat- 
ing from star formation processes (for a calculation using Nbody5 see [46], 
compare also [11,47] for the study of mass segregation in young forming star 
clusters by means of direct A^-body models). Finite sizes of stars lead to merg- 
ing in high-density phases and cause population gradients and unusually high 
frequencies of exotic objects like blue stragglers and pulsars in the cores and 
haloes of globular clusters. Attempts to model all these processes in direct 
A'-body models, with as many ingredients and realistic features included as 
possible are under way [2]. Ultimately we will be able from such models to 
provide synthetic observational data as e.g. color-magnitude diagrams. 
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